/*****************************************************************************

                                INIT.C

  Title:        General initialization package for SPEDEN's Solve.
  Author:       Hanna Szoke
  Date:         1/23/92
  Function:     This package contains set-up functions for Solve.
		These include: reading and interpreting .inp file,
		echoing what is in the input to the log file; and
                defining the run conditions for Solve.

		The basic .inp read function (used by ALL Speden's 
		pieces) is packaged in basics.c;
		Input of (hkl) data is in hklread.c.

******************************************************************************/

#include "speden.h"	/* includes util.h, cellparams.h, symmetry.h, dims.h */
#include "ccg.h"		/* needed for dfdxpc */

  /************************************************************
  declarations of some EXTERN variables in dims.h               
  ************************************************************/

real	nusqlim[NUMSEGS] ;	/* values for 1/8, 1/4,... of (hkl) shells */
real	fscale ;		/* scaling factor for fo data */
real	hrCutoff ;		/* hrCutoff * <vox> = high res limit	*/

  /************************************************************
  declarations of EXTERN variables in ccg.h               
  ************************************************************/

real	dfdxpc ;		/* frac. change, df/dx passed to getsol */

  /************************************************************
  declarations of EXTERN variables in speden.h               
  ************************************************************/

float	R_stop ;		/* R_factor at which code should stop	*/
real  discrp_frac ;           /* fraction of DISCRP set by user (def: 1) */

char    fo_filename[MAXSTRING] ;/* name of file, observed structure factors */

char	md_filename[MAXSTRING] ;/* name of file, model electron map */

int	Nconstraints ;		/* number of all constraints */
int	Ntargets ;		/* number of target constraints */
int	con_type[MAXCONSTR] ;	/* constraint type */
real	relwt_con[MAXCONSTR] ;	/* relative weight for a constraint */
real	cost_addend[MAXCONSTR] ;/* constant factor in the cost function */
real	con_coef[MAXCONSTR] ;	/* coefficient for cost function calc. */

char	ta_filename[MAXCONSTR][MAXSTRING] ;/* target file, phase_ext or target*/
char	wt_filename[MAXCONSTR][MAXSTRING] ;/* weight file, phase_ext or targe */
char	fc_filename[MAXCONSTR][MAXSTRING] ;/* fc list file, phase extension */
real	phase_ext_res ;		/* low resolution for phase extension */
int	thereis_phext ;
int	robust ;		/* if TRUE use robust hkl cost function */
real	rob_factor = UNSET ;	/* factor to use with robust cost function */

 /************************************************************
  declarations of EXTERN variables in util.h               
  ************************************************************/

char	sf_filename[MAXSTRING] ;	/* generic structure factor file name */
char	pwd_name[MAXSTRING] ;	/* name of str. factor file, w/o path */

  /************************************************************
  declarations of miscellaneous variables local to init.c  
  ************************************************************/

static	char	title[MAXSTRING] ;	/* whatever the user wants ... */

static	void	setGlobalDefaults() ;
static	void	readConstraintInfo() ;

				/*********************************************
               			Read Solve input 
				*********************************************/
void readInput()    
{
     int	k ;

     /************************************************************
     Read basics; set defaults for global variables.
     ************************************************************/

     readBasicInput() ;
     setGlobalDefaults() ;

     for (k = 0; k < Nlines; k++) {
          strcpy(nextline, allinp + k*MAXSTRING) ;
          sscanf(nextline, "%s", id) ;

          if (strcmp(id, "TITLE") == 0) {
               strcpy(title, nextline+6) ;
	       goodline[k] = TRUE ;
          }
          else if (strcmp(id, "DFDX_CRIT") == 0) {
               sscanf(nextline, "%*s %g", &t1) ;
               dfdxpc = t1 ;
	       goodline[k] = TRUE ;
          }

          else if ((strcmp(id, "FO_FILENAME") == 0) ||
              (strcmp(id, "FO_FN") == 0)) {
               sscanf(nextline, "%*s %s", fo_filename) ;
               goodline[k] = TRUE ;
          }
          else if (strcmp(id, "R_STOP") == 0) {
               sscanf(nextline, "%*s %g", &t1) ;
	       R_stop = t1 ;
               goodline[k] = TRUE ;
          }
          else if (strcmp(id, "DISCRP_FRAC") == 0) {
               sscanf(nextline, "%*s %g", &t1) ;
               discrp_frac = t1 ;
               goodline[k] = TRUE ;
   	  }
          else if (strcmp(id, "FSCALE") == 0) {
               sscanf(nextline, "%*s %g", &t1) ;
               fscale = t1 ;
               goodline[k] = TRUE ;
          }
          else if (strcmp(id, "ROB_FACTOR") == 0) {
               sscanf(nextline, "%*s %g", &t1) ;
               rob_factor = t1 ;
               goodline[k] = TRUE ;
          }
          if ((strcmp(id, "MD_FILENAME") == 0) ||
              (strcmp(id, "MD_FN")       == 0)) {
               sscanf(nextline, "%*s %s", md_filename) ;
               goodline[k] = TRUE ;
          }
          else if (strcmp(id, "NCONSTRAINTS") == 0) {
               sscanf(nextline, "%*s %d", &Nconstraints) ;
               goodline[k] = TRUE ;
          }
          else if ((HighRes) && (strcmp(id, "HRCUTOFF") == 0)) {
               sscanf(nextline, "%*s %g", &t1) ;
               hrCutoff = t1 ;
               goodline[k] = TRUE ;
          }
     }

     robust = fetchTFStatus("ROBUST", FALSE) ;
     if ((robust) && (rob_factor == UNSET))
         Fatal("Missing rob_factor!") ;

     if (strcmp(fo_filename, "none") == 0) 
         Fatal("Missing fo file name!") ;
          
     if (fscale == 0.)
         Fatal 
	     ("Missing fscale, scale factor for fo; when in doubt, use 1.") ;

     if (strcmp(md_filename, "none") == 0) 
	     Fatal("The starting binary file name is missing!");


     /************************************************************
     Deal with constraints.
     ************************************************************/

     readConstraintInfo() ;

     e_free(MAXSTRING*sizeof(char), allinp) ; 
}

				/*********************************************
               			Read Solve constraint information
				*********************************************/
void readConstraintInfo()    
{
     char	as_con_type[MAXSTRING][MAXCONSTR] ;
     int	j, k, m ;
     int	status = 0 ;

     for (k = 0; k < MAXCONSTR; k++) {
          strcpy(wt_filename[k], "none") ;
          strcpy(ta_filename[k], "none") ;
          strcpy(fc_filename[k], "none") ;
          strcpy(as_con_type[k], "") ;
	  relwt_con[k] = 0 ;
	  cost_addend[k] = 0 ;
	  con_coef[k] = 1 ;
     }
     thereis_phext = FALSE ;
     phase_ext_res = 0 ;

     if (Nconstraints == 0)
	  return ;

     for (k = 0; k < Nlines; k++) {
          strcpy(nextline, allinp + k*MAXSTRING) ;
          sscanf(nextline, "%s", id) ;

	/* General constraint input */

          if ((strncmp(id, "RELWT_CON", 9) == 0) && isdigit(id[9]))  {
	       j = atoi(&id[9]) - 1 ;
               sscanf(nextline, "%*s %g", &t1) ;
	       relwt_con[j] = t1 ;
	       goodline[k] = TRUE ;
          }
          else if ((strncmp(id, "COST_ADDEND", 11) == 0) && isdigit(id[11]))  {
	       j = atoi(&id[11]) - 1 ;
               sscanf(nextline, "%*s %g", &t1) ;
	       cost_addend[j] = t1 ;
	       goodline[k] = TRUE ;
	  }
          else if ((strncmp(id, "CON_TYPE", 8) == 0) && isdigit(id[8]))  {
	       j = atoi(&id[8]) - 1 ;
               sscanf(nextline, "%*s %s", as_con_type[j]) ;
	       goodline[k] = TRUE ;
	  }

	/* Target and Phase Extension constraint input */

          else if ((strncmp(id, "TA_FILENAME", 11) == 0) && isdigit(id[11])) {
	       j = atoi(&id[11]) - 1 ;
               sscanf(nextline, "%*s %s", ta_filename[j]) ;
	       goodline[k] = TRUE ;
	  }
          else if ((strncmp(id, "TA_FN", 5) == 0) && isdigit(id[ 5])) {
	       j = atoi(&id[5]) - 1 ;
               sscanf(nextline, "%*s %s", ta_filename[j]) ;
	       goodline[k] = TRUE ;
	  }
          else if ((strncmp(id, "WT_FILENAME", 11) == 0) && isdigit(id[11])) {
	       j = atoi(&id[11]) - 1 ;
               sscanf(nextline, "%*s %s", wt_filename[j]) ;
	       goodline[k] = TRUE ;
	  }
          else if ((strncmp(id, "WT_FN", 5) == 0) && isdigit(id[ 5])) {
	       j = atoi(&id[5]) - 1 ;
               sscanf(nextline, "%*s %s", wt_filename[j]) ;
	       goodline[k] = TRUE ;
	  }
          else if ((strncmp(id, "FC_FILENAME", 11) == 0) && isdigit(id[11])) {
	       j = atoi(&id[11]) - 1 ;
               sscanf(nextline, "%*s %s", fc_filename[j]) ;
	       goodline[k] = TRUE ;
	  }
          else if ((strncmp(id, "FC_FN", 5) == 0) && isdigit(id[ 5])) {
	       j = atoi(&id[5]) - 1 ;
               sscanf(nextline, "%*s %s", fc_filename[j]) ;
	       goodline[k] = TRUE ;
	  }

	/* Phase Extension resolution input */

          else if (strcmp(id, "PHASE_EXT_RES") == 0) {
               sscanf(nextline, "%*s %g", &t1) ;
	       phase_ext_res = t1 ;
	       goodline[k] = TRUE ;
          }
     }

     /************************************************************
     Check each constraint type for legality.
     ************************************************************/

     for (k = 0; k < Nconstraints; k++) 
          if (relwt_con[k] == 0) {
               sprintf(message, "Missing relative weight for target # %d!", k+1) ;
	       Warning(message) ;
               status = -1 ; 
           }

     for (k = 0; k < Nconstraints; k++) {
	  
	  for (m=0; m<(int)strlen(as_con_type[k]); m++)
	       as_con_type[k][m] = tolower((int) as_con_type[k][m]) ;

	  if (strcmp(as_con_type[k], "target") == 0) {
	       con_type[k] = TARGET ;
	       Ntargets++ ;
               if (strcmp(ta_filename[k], "none") == 0) {
	           sprintf(message, "Missing ta_filename for target # %d!", 
		   k+1) ;
	           Warning(message) ;
	           status = -1 ; 
               }
          }
	  else if (strcmp(as_con_type[k], "phase_ext") == 0) {
               thereis_phext = TRUE ;
	       con_type[k] = PHASE_EXT ;
               ph_coef = relwt_con[k] ;
               if (strcmp(fc_filename[k], "none") == 0) {
	           sprintf(message, 
	           "Missing filename in reciprocal space for phase extension!") ; 
	           Warning(message) ;
	           status = -1 ; 
               }
	       if (phase_ext_res == 0) {
		   Warning("Missing phase extension resolution!") ;
	           status = -1 ; 
               }
	       else if (phase_ext_res < input_res) {
		   sprintf(message,
	           "Phase_ext_res (%g A) should be no higher than input_res (%gA)!",
	           phase_ext_res, input_res) ;
	           Warning(message) ;
	           status = -1 ; 
               }
          }
	  else if (strcmp(as_con_type[k], "cs") == 0) {
	       con_type[k] = CS ;
	       ;	/* no special checks */
          }
	  else if (strcmp(as_con_type[k], "sayre") == 0) {
	       Warning( 
	  "There is no Sayre's equation cost function in Speden!") ;
	       status = -1 ; 
          }
	  else  {
	       Warning( "Illegal constraint type") ;	
	       status = -1 ; 
          }
     }
     if (status != 0)
          Fatal("1 or more illegal constraint entries.") ;
}


#define	DFDX_SOLVE	3.e-2 
				/*********************************************
                       		Set defaults for global variables  
				*********************************************/
void setGlobalDefaults()    
{
     dfdxpc = DFDX_SOLVE ;		

     R_stop = 0. ;
     discrp_frac = 1. ;
     fscale = 0. ;
     Nconstraints = 0 ;
     Ntargets = 0 ;
     hrCutoff = 10. ;

     strcpy(fo_filename, "none") ;
     strcpy(md_filename, "none") ;
     strcpy(title, "") ;
}

				/*********************************************
          			Echo information in .inp file to log.     	
				*********************************************/
void echo()      
{
     int	k ;

     if ((int)strlen(title) > 0) 
          fprintf(fp_log, "\n%s\n\n", title) ;
     
     /************************************************************
     Information regarding input files.
     ************************************************************/

     fprintf(fp_log, "\n") ;
	  if ((int)strlen(fo_filename) < 40)
               fprintf(fp_log, 
        "Observed structure factors will be read from %s\n", fo_filename) ;
	  else {
               fprintf(fp_log, 
	       "Observed structure factors will be read from\n") ;
               fprintf(fp_log, " %s\n", fo_filename) ;
          }
          if (useSig) 
	       fprintf(fp_log, "Sigmas will be used for weighting\n") ;
          else
	       fprintf(fp_log, "Sigmas will not be used for weighting\n") ;
	  fprintf(fp_log, "Data scaling factor is %g\n", fscale) ;

     if ((int)strlen(md_filename) < 40)
          fprintf(fp_log, 
     "Starting physical space model will be read from %s\n", md_filename) ;
     else 
          fprintf(fp_log, 
     "Starting physical space model will be read from\n %s\n", md_filename) ;
     
     /************************************************************
     Constraints, relative weights except for constraints already
     reported in the mode-of-operation banner. 
     ************************************************************/

     fprintf(fp_log, "\n") ;
     fprintf(fp_log, "Relative weight for reciprocal space is 1. \n") ;

     fprintf(fp_log, "\n") ;

     for (k = 0; k < Nconstraints; k++) {
	  switch (con_type[k]) {
		  case TARGET: 
	       if ((int)strlen(ta_filename[k]) < 40)
                    fprintf(fp_log, 
                  "Target Np values will be read from %s\n", ta_filename[k]) ;
	       else {
                    fprintf(fp_log, "Target Np values will be read from\n") ;
                    fprintf(fp_log, " %s\n", ta_filename[k]) ;
               }

	       if ((int)strlen(wt_filename[k]) < 40)
                    fprintf(fp_log, 
                  "Target weights will be read from %s\n", wt_filename[k]) ;
	       else {
                    fprintf(fp_log, "Target weights will be read from\n") ;
                    fprintf(fp_log, " %s\n", wt_filename[k]) ;
               }

               fprintf(fp_log, 
	         "Target relative weight is %g \n", relwt_con[k]) ;
               break ;
     
	  default :
               break ;

          }
     }

     fprintf(fp_log, "\n") ;
     fflush(fp_log) ;
}

void	echo_ignored()
{
     FILE	*fp ;
     int	k = 0 ;
     char 	first = TRUE ;

     /* Was anything been ignored? */

     if ((fp = fopen(input_filename, "r")) == NULL)
     {
          sprintf(message, "Cannot open %s", input_filename);
          Fatal(message) ;
     }

     while (fgets(nextline, MAXSTRING, fp) != NULL)  {
	  if (net_strlen(nextline) > 1) {
	       if (!goodline[k]) {
                    if (first) {
                         printTwice(
                    "The following input lines were ignored:\n") ;
                         first = FALSE ;
                    }
	            sprintf(message, "%s%s", ">>  ", nextline) ;
	            printTwice(message) ;
               }
               k++ ;
          }
     }
     printTwice("") ;
     fclose(fp) ;
}
